RStudio is a GUI version of R, which is command-line only. RStudio includes the following:

• Manual, help function
• Script separate from command-line
• Lists your variables
• Easy install of new packages
• Plots are shown within RStudio

If you work with linux (or want to work from with R from the command line but still document), start R by opening a terminal, changing directory to the R_exercises folder and then just typing R. Then, you should see something like this:

Below is how to run it in a terminal (i.e. on ada). For this to work just remove the # symbol.

#Ask what R version we have
#R.version

#start R
#R

#exit R
#q()

General

One neat way to work with R is to combine R-code with text in R markdown format (rmd). The R markdown file specifies code chunks which will be executed in R (or python or bash) and plain text which will be written to the report as is. The report is created by rendering the file in R, then the R-code is executed and the results are merged in the pdf or html output.

The easiest way to get started is to create a first R markdown document using Rstudio: Go to the menu and press “File/new File/R markdown”. This will create an R markdown file that already contains some R code and text. Do this, and look at the contents of this document.

Two useful information sources are accessible from the Rstudio menu: help/cheatsheets/R markdown Cheat Sheet and help/cheatsheets/R markdown Reference Guide. Short help on the text format that can be used in markdown documents can be found in help/Markdown quick reference.

If you work with plain R you can open a new file in any text editor and save it with the .rmd extension and research the format requirements on your own. Check out the rmarkdown website and the cheat sheet in the R_exercises folder to get started on that (and read on). You can find more R cheat sheets on https://rstudio.com/resources/cheatsheets/

R code in R markdown

• The R-code is embedded inbetween the “‘{r} and “‘ symbols. In Rstudio, on the top-right position of such a section you will find three symbols. Pressing the middle one will run all code chunks above, while the right symbol will run the current R-chunk.
• An important menu-button is the “Knit” button. Pressing this button will create the ultimate document. You can choose between HTML, PDF and WORD documents.

The cheat sheet can be found here: https://www.rstudio.com/wp-content/uploads/2016/03/rmarkdown-cheatsheet-2.0.pdf and will help you also if you want to modify the treatment of code chunks. For example, per default, the R code and its output are written to the report, but you can suppress one or both of that.

Text in markdown

• Headings are defined with “#”, or “##”, or “###” for first, second and third level.
• Lists are created by using “*” for (bullets) and “1”, … for numbered lists.
• A lot of other options are possible, see the cheat sheets.

Markdown options

• include = FALSE — prevents code and results from appearing in the finished file. R Markdown still runs the code in the chunk, and the results can be used by other chunks.
• echo = FALSE —prevents code, but not the results from appearing in the finished file. This is a useful way to embed figures.
• message = FALSE —prevents messages that are generated by code from appearing in the finished file.
• warning = FALSE —prevents warnings that are generated by code from appearing in the finished.
• fig.cap = “…” —adds a caption to graphical results.

1. Online

• www.r-project.org stack overflow many more

2. Inside of R

help(mean)
?mean

The workspace is your current R working environment and includes any user-defined objects (vectors, matrices, data frames, lists, functions). At the end of an R session, the user can save an image of the current workspace that is automatically reloaded the next time R is started.

• From now on the code is in the grey area.
• If you see a line symbol, such as [1] then this shows the result of running our code. For example, first we run the code getwd() and the second box shows us the results.

#print the current working directory
getwd() 

[1] "/Users/ninadombrowski/Desktop/WorkingDir/Github_Workflows/R_tutorial"

#list the objects in the current workspace
ls()   

[1] "encoding"      "growth_data"   "input_file"    "out_dir"      
[5] "test_variable"

A unit that you work with, i.e. data and functions.

Characteristics = Data type (kind of data) and data structure (the container)

• Data type = numeric, characters, factors, … You can identify your data type with the function typeof()

• Data structure= vector, matrix, list. You can identify the structure with class() or is()

• str() = defines what the indiv. components of your data are

The directory from which you work is usually first set from where you start R. But it can be re-set to find your data more easily. Ideally, you make one wdir per project and define the path in the script (see later below). It is recommended to have a similar format for these project folders, i.e. consider to create subfolders for input and output files. From the wdir you set you can load files using absolute and relative paths.

An example would be something with a structure like this:

Here, we have 4 projects, and in each folder we have the R_script and folders for the required input and required output files. Also useful to have is a text file with the session info and if you plan to share this with others it is also good to have a README file that provides some background on the analysis done.

Options to see the working dir and set the working directory in R are:

getwd()

setwd(getwd())

Packages are a collection of functions and tools that can be added to R that are often contributed by the community. Careful: there might be incompatibilites and packages are updated frequently but updating can break dependencies.

You need to install packages and load them EVERY TIME you want to use them. Therefore, ideally add them at the beginning of your scripts.

Installing packages

We have two ways to install packages:

1. Via the console by typing:

install.packages(“package-name”)

This will download a package from one of the CRAN mirrors assuming that a binary is available for your operating system. If you have not set a preferred CRAN mirror in your options(), then a menu will pop up asking you to choose a location.

Use old.packages() to list all your locally installed packages that are now out of date. update.packages() will update all packages in the known libraries interactively. This can take a while if you haven’t done it recently. To update everything without any user intervention, use the ask = FALSE argument.

2. Using R studio:

Go to the lower right hand-side window and click on packages and then install. Find the packages your are interested in.

Notice: If libraries come with their own data (i.e. example tables), then the data needs to be loaded separately. I.e. via data(cars)

Loading packages into your current R session

R will not remember what libraries you have loaded after you closed R. Therefore, you need to load necessary libraries everytime. Here, we will load some libraries that are usually quite helpful and it is recommended to make the libraries you load part of each of your scripts.

#needed for Rmarkdown
library(knitr)
library(kableExtra)

#other useful packages for data treatment, plotting and stats
library(tinytex)
library(stats)
library(knitr)
library(kableExtra)
library(ggplot2)
library(data.table)
library(tidyr)
library(plyr)

The <- symbol assigns a value to a variable,

General rules for the syntax:

• If you work with characters, i.e. words like ‘hello’, then this need to written like this: “hello” (this will become clearer below)
• R is case sensistive
• If a variable exists, it will overwrite it with a new variable without asking
• ls() shows all the variables that are known by the system at the moment
• you can remove indiv. objects with rm() and remove all objects rm(list=ls())

For example:

We can store our path for the working directory in a variable, which is useful if we want to re-run our script from the same folder.

Notice can the path to whereever you download the tutorial.

#store a path (a character string) to a variable with the name wdir
wdir <- "~/Desktop/WorkingDir/Github_Workflows/R_tutorial"

#print variable to screen
wdir

[1] "~/Desktop/WorkingDir/Github_Workflows/R_tutorial"

#use our new  variable to set our wdir
setwd(wdir)

We can store more or less everything in a variable and use it later.

#store some numbers
x <- 1
y <-4

#do some simple math with the numbers we have stored
x+y

[1] 5

Functions are build in code that we can use to make our life easier, i.e. we can calculate lengtsh of vectors or do math and do statistical analyses. Base R already knows many useful functions but loading new packages greatly increases our repertoire. A list of most used functions can be found here: http://www.statmethods.net/management/functions.html

A function consists of:

1. Function name
2. Arguments (optional, some might be set with a default)
3. Body of the function = defines what the function does
4. return = in case the function returns a variable

As an example lets test some simple functions: print and log:

#use the print function
print(3+5)

[1] 8

#use the log function
log(10)

[1] 2.302585

Call the default values of a function

Every function commes with a set of arguments that you can set but that usually also have some default values. In R Studio you can easily access all those details with the help function.

• ? allows us to first of all check exactly what a function is doing. If you scroll down to the bottom of the help page you also get some examples on how to use a function.
• More specifically the help function also allows us to get details on the arguments of a function.
• For example, if we check the help page of read.table we see that by default this function does not read in a header and if we want to provide a header we have to change that argument.

#let's check what **log** is doing
?log

#lets check the default arguments of a function
?read.table

Other useful functions:

To work with any kind of data we need to first read the data into R to be able to work with it.

For tables, there are some things to be aware of:

• It matters what separator our table uses to specificy individual columns. I.e. some programs store data using commas while others use tab as default delimiter. By default R assume we use a tab, but we can change this behaviour when we read in the table.
• Do not have any hash symbols (#) in your table. R will read this as a commented cell and not read the table from that point onward
• Avoid empty cells, as these sometimes can mess up your data.

For now, let’s read in the table with our growth data and store it under the variable name growth_data

Notice, to read in the file GenomeInfo.txt we need to direct it to the correct path as we do not have the file in the working directory but in a subdirectory.

Options that are good to leep in mind when reading in a table:

• sep = define that our field separator is a tab. A tab is written like this /t. If your data is using a space or comma, you can change that here.
• header = tell R that our data comes with a custom header (the first row in your datafra)
• quote = deals with some annoying issue with data formatting in excel files

General notice:

• To view data, the head() command is extremely practical, use it always when you modify data to check if everything went alright
• dim() is another useful function that displays the dimensions of a table, i.e. how many rows and columns we have. Again, this is useful to verify our data after we have transformed it to check if everything went alright.
• colnames() allows to only we the column names
• rownames() allows to only we the row names. Usually these are numbers, but we can also add anything else into the rows.

#read in data
growth_data <- read.table("Input_files/Growth_Data.txt", sep="\t", header=T,  quote = "")

#check the first few lines of our data
kable(head(growth_data), format='markdown')

#check the dimensions of our data
dim(growth_data)

[1] 105   5

#check the column names
colnames(growth_data)

[1] "SampleID"      "Nutrient"      "Condition"     "FW_shoot_mg"  
[5] "Rootlength_cm"

#check the row names
rownames(growth_data)

  [1] "1"   "2"   "3"   "4"   "5"   "6"   "7"   "8"   "9"   "10"  "11"  "12" 
 [13] "13"  "14"  "15"  "16"  "17"  "18"  "19"  "20"  "21"  "22"  "23"  "24" 
 [25] "25"  "26"  "27"  "28"  "29"  "30"  "31"  "32"  "33"  "34"  "35"  "36" 
 [37] "37"  "38"  "39"  "40"  "41"  "42"  "43"  "44"  "45"  "46"  "47"  "48" 
 [49] "49"  "50"  "51"  "52"  "53"  "54"  "55"  "56"  "57"  "58"  "59"  "60" 
 [61] "61"  "62"  "63"  "64"  "65"  "66"  "67"  "68"  "69"  "70"  "71"  "72" 
 [73] "73"  "74"  "75"  "76"  "77"  "78"  "79"  "80"  "81"  "82"  "83"  "84" 
 [85] "85"  "86"  "87"  "88"  "89"  "90"  "91"  "92"  "93"  "94"  "95"  "96" 
 [97] "97"  "98"  "99"  "100" "101" "102" "103" "104" "105"

Useful comments:

Sometimes we have to deal with really large data that take long to load with read.table. The function fread() from the data.table package is a very nice alternative. Below, we will explain how to read in R libraries to use additional functions.

This script always uses kable for plotting, which is just to make it visually attractive in html. Whenever you see a function using kable simply replace it with the head() function.

Now, if we would have modified the table we might want to store it on your computer. We can do this using write.table() and below we use a different output directory. Notice, we always start from the location we set as working directory.

write.table(growth_data, "Output_examples/growth_data_changed.txt",  sep = "\t", row.names = T, quote =F)

• sep –> we define what delimiter we want to use
• row.names = T –> we want to include whatever is in the rownames
• quote = F –> we do not want R to add any quotes around our columns.

A vector is a collection of items of the same type (i.e characters, numbers). You can read in numbers and characters into the same vector, however, the number will be then seen as a character.

• length() = a function that allows us to quickly ask how long an object, such as a vector, is.
• c()= a function that will create a vector (a one dimensional array) and in our case store 3 numbers. We need to use this everytime we deal with more than one number, character, etc….
• square brackets [] = allow us to retrieve certain elements of a vector, i.e. [3] retrieves the 3rd element
• we can combine c() and [] if we want to retrieve several elements of a vector.

#lets create a random vector
our_own_vector <- c(2, 3, 5) 

#show the vector we just created
our_own_vector

[1] 2 3 5

#retrieve the third element stored in a vector
our_own_vector[3]

[1] 5

#retrieve the 1st and 3rd element by combining ``c()`` and []
our_own_vector[c(1,3)]

[1] 2 5

#asking how long your vector is
length(our_own_vector)

[1] 3

#we can also add vectors of the same length together
x <- c(1,2,3,4)
y <- c(1,2,3,4)

#and now we can combine our vectors
x + y

[1] 2 4 6 8

Beware: If we add two vectors of different length, the shorter vector is duplicated. This only works if the shorter vector is proportional to the longer one

#adding vectors of different lengths
x <- c(1,2)
y <- c(1,2,3,4)

#and now we can combine our vectors
x + y

[1] 2 4 4 6

Vector indexing

If we want to only retrieve part of the data stored in a vector we can create a subset using the index between square brackets [], e.g.,

#create a vector
x <- 1:10
x

 [1]  1  2  3  4  5  6  7  8  9 10

#retrieve the second element from our vector
x[2]

[1] 2

#retrieve the 4-6th element
x[4:6]

[1] 4 5 6

We not only can extract the nth element but if we have header names then we can also use these to retrieve data:

#create a vector and give it names
x <- c(3,4,5)
names(x) <- c("oranges","apples","bananas")

#check how our data looks
x

oranges  apples bananas 
      3       4       5 

#now we can retrieve part of the vector using the names
x[c("apples","oranges")]

 apples oranges 
      4       3 

We can also change elements in our vector:

#create a vector
x <- 1:10

#change the second last positions to 5 and 9
x[9:10] <- c(5,9)

#check if this worled
x

 [1] 1 2 3 4 5 6 7 8 5 9

#we can not only add things, we can also remove this using the minus/**-** symbol
x[-3]

[1] 1 2 4 5 6 7 8 5 9

#if we want to remove more than one thing we can use the **c()**
x[-(4:10)]

[1] 1 2 3

Dealing with NAs in our data

#create a vector that includes an NA
y <- c(1,2,3,NA,5)

#check whether we have NAs
is.na(y)

[1] FALSE FALSE FALSE  TRUE FALSE

#remove NAs from our data
y[!is.na(y)]

[1] 1 2 3 5

#create a vector that includes an NA
x <- c(1,2,3,NA,5)

#check whether we have NAs
is.na(x)

[1] FALSE FALSE FALSE  TRUE FALSE

#remove NAs from our data
x[is.na(x)] <- 0

#check data
x

[1] 1 2 3 0 5

Matrices are the R objects in which the elements are arranged in a two-dimensional rectangular layout. They contain elements of the same type. Although you can construct matrices with characters or logicals, matrices are generally used to store numeric data. The basic syntax for creating a matrix is:

matrix(data, nrow, ncol, byrow, dimnames)

• data: input vector whose components become the data elements from the matrix.
• dnrow: number of rows to be created.
• dncol: number of columns to be created.
• dbyrow: logical. If FALSE,(the default) the matrix is filled by columns, otherwise the matrix is filled by rows.
• ddimnames: A ‘dimnames’ attribute for the matrix: NULL or a list of length 2 giving the row and column names respectively.

In contrast in a data frame (see below) the columns contain different types of data, while in a matrix all the elements are the same type of data. A matrix in R is like a mathematical matrix, containing all the same type of thing (usually numbers). R often but not always can use dataframes and a matrix used interchangably.

• Individual elements in a matrix can be printed using [row,column]. For example [2,3] would pull out the value in the 2nd ROW and third COLUMN.
• dim() is extremely useful to control whether our data was transformed correctly during different operations. For example, after we merge two files we would like to know that they still have the same number of rows as when we started the analysis. Same if we remove for example 10 samples, then we want to make sure that this is indeed what happened.

#define our row and column names
row.names = c("row1", "row2", "row3", "row4")
col.names = c("col1", "col2", "col3")

#create our matrix (check the help function to see what is happening)
matrix_A <- matrix(c(3:14), nrow = 4, byrow = T, dimnames = list(row.names,col.names))

#check how our matrix looks like
matrix_A

     col1 col2 col3
row1    3    4    5
row2    6    7    8
row3    9   10   11
row4   12   13   14

#print the value in the 2row and 3rd column
matrix_A[2,3]

[1] 8

#print the values in the 3rd column
matrix_A[,3]

row1 row2 row3 row4 
   5    8   11   14 

#print everything except the 1st row
matrix_A[-1,]

     col1 col2 col3
row2    6    7    8
row3    9   10   11
row4   12   13   14

#print everything except the 2nd column
matrix_A[,-2]

     col1 col3
row1    3    5
row2    6    8
row3    9   11
row4   12   14

#see the dimensions of matrix, i.e. the nr of rows and columns
dim(matrix_A)

[1] 4 3

Sometimes you need to store data of different types. For example, if you are collecting data of a patient, you might want to have cell counts (numeric), the microbes we investigated (character), their status (logical, with TRUE for alive and FALSE for dead, for instance), etc. This kind of data can be stored in lists. Lists are the R objects which contain elements of different types (numeric, strings, vectors, even another list, or a matrix).

A list is created using the list() function.

For example, the following variable x is a list containing copies of three vectors n, s, b, and a numeric value 3 (which might be the replicate number of something else).

#define our vectors
n = c(20, 30, 50) 
s = c("Ecoli", "Archaeoglobus", "Bacillus") 
b = c(TRUE, FALSE, TRUE) 

#combine the vectors in a list
our_list = list(counts=n, strain=s, status=b, random = 3) 

#show our list
our_list

$counts
[1] 20 30 50

$strain
[1] "Ecoli"         "Archaeoglobus" "Bacillus"     

$status
[1]  TRUE FALSE  TRUE

$random
[1] 3

#sublist the second element in a list
our_list[2]

$strain
[1] "Ecoli"         "Archaeoglobus" "Bacillus"     

#retrieve the 2nd and 3rd member of our list
our_list[c(2, 3)] 

$strain
[1] "Ecoli"         "Archaeoglobus" "Bacillus"     

$status
[1]  TRUE FALSE  TRUE

#we can also retrieve elements of a list if we know the name using two different ways:
our_list$strain

[1] "Ecoli"         "Archaeoglobus" "Bacillus"     

our_list[["strain"]]

[1] "Ecoli"         "Archaeoglobus" "Bacillus"     

In the last example we use the dollar symbol to extract data, i.e. to extract variables in a dataset (a matrix, list, dataframe). I.e. above the data we want to access is ‘our_list’ and the variable we want to extract is the strain.

The R-Object data frame is a table in which each column contains values of one variable type and each row contains one set of values from each column. You can think of a data frame as a list of vectors of equal length. Most of our data very likely will be stored as dataframes.

A Dataframe usually follows these rules:

• Column names (i.e. the header of our data) should be non-empty (if they are, R provides the object with default values).

• Row names should be unique

• Each column should contain the same number of data items

• Numeric and logical vectors and factors are included as is, whereas character vectors are coerced to factors, unless one sets the argument stringsAsFactors to F.

• typeof(df.object) is list (type or storage mode) and class(df.object) is data.frame (the class of the object).

• The top line of the table, called the header, contains the column names.

• Each horizontal line afterward denotes a data row, which begins with the name of the row, and then followed by the actual data.

• Each data member of a row is called a cell.

Importantly, most of the things we have learned before, i.e. how to subset data, apply here too.

The growth data we provide as test data is a dataframe. Therefore, we can check this data for different characteristics of a dataframe.

Viewing data Dataframes

• We can use the brackets as before to extract certain rows or columns.
• We can use the dollar sign to again extract information as long as we know the column names. I.e. now we want to access the shoot fresh weight (FW_shoot_mg) in our ‘growth_data’ dataframe.
• colnames() allows us to access the column names, i.e. the headers
• rownames() allows us to access the rownames of our data (usually these are numbered if not specified otherwise while reading the table)
• dim() allows us to check the dimensions (i.e. the number of rows and columns). This is useful to regullary check, especially if we modified our data somehow.

#view our table
kable(head(growth_data), format='markdown')

#extract the data from the 2nd row
growth_data[2,]

  SampleID Nutrient Condition FW_shoot_mg Rootlength_cm
2      noP      noP      MgCl        6.52      5.743447

#extract the first three columns
head(growth_data[,1:3])

  SampleID Nutrient Condition
1      noP      noP      MgCl
2      noP      noP      MgCl
3      noP      noP      MgCl
4      noP      noP      MgCl
5      noP      noP      MgCl
6      noP      noP      MgCl

#extract a column of our data using the column name
head(growth_data$FW_shoot_mg)

[1] 10.26  6.52 12.17 11.37  9.80  3.75

#print our headers
colnames(growth_data)

[1] "SampleID"      "Nutrient"      "Condition"     "FW_shoot_mg"  
[5] "Rootlength_cm"

#print the rownames
rownames(growth_data)

  [1] "1"   "2"   "3"   "4"   "5"   "6"   "7"   "8"   "9"   "10"  "11"  "12" 
 [13] "13"  "14"  "15"  "16"  "17"  "18"  "19"  "20"  "21"  "22"  "23"  "24" 
 [25] "25"  "26"  "27"  "28"  "29"  "30"  "31"  "32"  "33"  "34"  "35"  "36" 
 [37] "37"  "38"  "39"  "40"  "41"  "42"  "43"  "44"  "45"  "46"  "47"  "48" 
 [49] "49"  "50"  "51"  "52"  "53"  "54"  "55"  "56"  "57"  "58"  "59"  "60" 
 [61] "61"  "62"  "63"  "64"  "65"  "66"  "67"  "68"  "69"  "70"  "71"  "72" 
 [73] "73"  "74"  "75"  "76"  "77"  "78"  "79"  "80"  "81"  "82"  "83"  "84" 
 [85] "85"  "86"  "87"  "88"  "89"  "90"  "91"  "92"  "93"  "94"  "95"  "96" 
 [97] "97"  "98"  "99"  "100" "101" "102" "103" "104" "105"

#get the dimensions
dim(growth_data)

[1] 105   5

Adding new columns to Dataframes

Below is a very basic way to add a new column (we name it newColumn) and fill all rows with the word Experiment1

#expand a dataframe, functions data.frame or cbind (or see below)
growth_data$newColumn <- "Experiment1"

There are more sophisticated ways to add columns based on conditions or even merge dataframes. All of these things we will discuss later.

If we have our own table, we might want to see as what kind of data our table is stored. We have several ways to do this

• class() = determines as what kind of object is stored
• str() = display the internal structure of an R object.

There are different structures that can be stored in a dataframe.

1. character: “a”, “swc”
2. numeric: 2, 15.5
3. integer: 2L (the L tells R to store this as an integer)
4. logical: TRUE, FALSE
5. complex: 1+4i (complex numbers with real and imaginary parts)

#check what kind of data we have:
class(growth_data)

[1] "data.frame"

#check how are different parts of our data stored?
str(growth_data)

'data.frame':   105 obs. of  6 variables:
 $ SampleID     : Factor w/ 8 levels "noP","noP_101",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Nutrient     : Factor w/ 2 levels "noP","P": 1 1 1 1 1 1 1 1 1 1 ...
 $ Condition    : Factor w/ 4 levels "MgCl","Strain101",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ FW_shoot_mg  : num  10.26 6.52 12.17 11.37 9.8 ...
 $ Rootlength_cm: num  5.93 5.74 6.83 6.74 6.74 ...
 $ newColumn    : chr  "Experiment1" "Experiment1" "Experiment1" "Experiment1" ...

Factors are the data objects which are used to categorize the data and store it as levels. They can store both strings and integers. They are useful in the columns which have a limited number of unique values, e.g., “Male”, “Female”, TRUE, FALSE etc. They are also useful in data analysis for statistical modeling.

gender <- c("male","male","female", "female", "male","female")
gender

[1] "male"   "male"   "female" "female" "male"   "female"

#what class do we have
class(gender)

[1] "character"

typeof(gender)

[1] "character"

Here, we see that the variable gender is of class factor.

A level is a vector of all possible values of a certain column.

Check the structure of part of your dataframe with:

growth_data$Nutrient

  [1] noP noP noP noP noP noP noP noP noP noP noP noP noP noP noP noP noP noP
 [19] noP noP noP noP noP noP noP noP noP noP noP noP noP noP noP noP noP noP
 [37] noP noP noP noP noP noP noP noP noP noP noP noP noP noP noP noP P   P  
 [55] P   P   P   P   P   P   P   P   P   P   P   P   P   P   P   P   P   P  
 [73] P   P   P   P   P   P   P   P   P   P   P   P   P   P   P   P   P   P  
 [91] P   P   P   P   P   P   P   P   P   P   P   P   P   P   P  
Levels: noP P

Here, we can see that Nutrient data column has 105 elements and 2 levels (P and noP).

We can check the class of an indiv. variable using class() or of all variables of a dataframe using str() .

class(growth_data$Nutrient)

[1] "factor"

str(growth_data)

'data.frame':   105 obs. of  6 variables:
 $ SampleID     : Factor w/ 8 levels "noP","noP_101",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Nutrient     : Factor w/ 2 levels "noP","P": 1 1 1 1 1 1 1 1 1 1 ...
 $ Condition    : Factor w/ 4 levels "MgCl","Strain101",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ FW_shoot_mg  : num  10.26 6.52 12.17 11.37 9.8 ...
 $ Rootlength_cm: num  5.93 5.74 6.83 6.74 6.74 ...
 $ newColumn    : chr  "Experiment1" "Experiment1" "Experiment1" "Experiment1" ...

This prints our session info (used R version, packages, etc) and is very important to store to reproduce our code for others

sessionInfo()

R version 3.5.0 (2018-04-23)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13.6

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] plyr_1.8.6        tidyr_1.0.2       data.table_1.12.8 ggplot2_3.3.2    
[5] tinytex_0.24      kableExtra_1.1.0  knitr_1.29       

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.4        pillar_1.4.4      compiler_3.5.0    highr_0.8        
 [5] rmdformats_0.3.7  tools_3.5.0       digest_0.6.27     evaluate_0.14    
 [9] lifecycle_0.2.0   tibble_3.0.4      gtable_0.3.0      viridisLite_0.3.0
[13] pkgconfig_2.0.3   rlang_0.4.8       rstudioapi_0.11   yaml_2.2.1       
[17] xfun_0.15         withr_2.3.0       dplyr_1.0.2       stringr_1.4.0    
[21] httr_1.4.1        xml2_1.2.5        generics_0.0.2    vctrs_0.3.4      
[25] hms_0.5.3         tidyselect_1.1.0  webshot_0.5.2     grid_3.5.0       
[29] glue_1.4.2        R6_2.5.0          rmarkdown_2.3     bookdown_0.20    
[33] purrr_0.3.3       readr_1.3.1       magrittr_1.5      scales_1.1.1     
[37] ellipsis_0.3.0    htmltools_0.4.0   rvest_0.3.5       colorspace_1.4-1 
[41] stringi_1.5.3     munsell_0.5.0     crayon_1.3.4     

Notice: At the moment, we do not want to clean up the working directory, therefore, in our documentation we put a # in front of it. If you want to run the code just remove the #

#rm(list = ls()) 

One of the most important features of R (and of any programming language) are logical operators. When R evaluates an expression containing logical operators, it will return either TRUE or FALSE. It follows a list of them:

#basic operators
1==1

[1] TRUE

2>1

[1] TRUE

1==1 & 1>2

[1] FALSE

(1==1& 1>2) | 5<6

[1] TRUE

#only print data for the normal P conditions
P_data <- growth_data[growth_data$Nutrient == "P", ]

#test that all went ok
dim(growth_data)

[1] 105   6

dim(P_data)

[1] 53  6

#we can also subset our data by numbers, i.e. only keep rows were roots are longer than 10cm
growth_data_10cm <- growth_data[growth_data$Rootlength_cm > 10, ]
dim(growth_data_10cm)

[1] 47  6

#instead of keeping rows with roots longer than 10cm we can also remove them using *!*
growth_data_no10cm <- growth_data[!growth_data$Rootlength_cm > 10, ]
dim(growth_data_no10cm)

[1] 58  6

Above you see a good example why dim() can very useful as it allows us to quickly check that what we do actually works. Always check afterwards that your function did what you expect you to do since as long as R is able to do something it will do it regardless if that might be different that you want R to do.

If you have a file that stores R commands, e.g. custom functions or (parts of) an analysis, let’s say my_script.R in your working directory, they can be executed in the R session by typing the command source. For example, collaborators might give you scripts or you later might decide to write your own.

#source("my.script")

This section will be extended at a later point.

Apart from the elementary operations, all of the common arithmetic functions are available: log, exp, sin, cos, tan, sqrt, etc. Other useful functions one can use on vectors are:

Data that can be changed is called mutable, while data that cannot be is called immutable. Like strings, numbers are immutable: there’s no way to make the number 0 have the value 1 or vice versa (at least, not in R—there actually are languages that will let people do this, with predictably confusing results). Vectors, data frames, and matrices, on the other hand, are mutable: they can be modified after they have been created.

Programs that modify data in place can be harder to understand than ones that don’t because readers may have to mentally sum up many lines of code in order to figure out what the value of something actually is.

On the other hand, programs that modify data in place instead of creating copies that are almost identical to the original every time they want to make a small change are much more efficient. There are many ways to change the contents of a vector and some examples are listed below.

We also use a new function:

• append() –> Add elements to a vector.

#lets make a vector
our_own_vector <- c(1, 3, 5, 7, 9)
our_own_vector

[1] 1 3 5 7 9

#add another datapoint to our vector
our_own_vector <- append(our_own_vector, 13)
our_own_vector

[1]  1  3  5  7  9 13

#add +1 to all our four numbers
our_own_vector <- our_own_vector + 1
our_own_vector

[1]  2  4  6  8 10 14

#remove the first element of our vector
our_own_vector <- our_own_vector[-1]
our_own_vector

[1]  4  6  8 10 14

We already have seen important ways to subset data:

1. Use of the index by using the square brackets

#subsetting rows and columns using the index
growth_data[1:3,2:4]

  Nutrient Condition FW_shoot_mg
1      noP      MgCl       10.26
2      noP      MgCl        6.52
3      noP      MgCl       12.17

2. Use of operators

#subset using operators (only print rows if the Nutrient column equals P)
P_data <- growth_data[growth_data$Nutrient == "P", ]

grep and grepl

Oftentimes you may need to filter a data set based on a partial character string that is beyond the scope of comparison operators.

R provides such functions (grep and grepl) that match character patterns in specified vector. While both of these functions find patterns, they return different output types based on those patterns.

• grep returns numeric values that correspond to the indexed locations of the patterns
• grepl returns a logical vector in which “TRUE” represents a pattern match.

In our growth data, we only want to print measurements of our strains.

#grep return the index value of each matched pattern
grep("Strain", growth_data$Condition)

 [1]  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32
[20]  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51
[39]  52  68  69  70  71  72  73  74  75  76  77  78  79  80  81  82  83  84  85
[58]  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 101 102 103 104
[77] 105

#grepl returns a logical output for each indices in the original vector 
#with "TRUE" representing matched patterns
grepl("Strain", growth_data$Condition)

  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [13] FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
 [25]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
 [37]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
 [49]  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE
 [73]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
 [85]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
 [97]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

Now lets use this to actually filter our data table for a pattern.

filter_for_value <-growth_data[grepl("Strain", growth_data$Condition),]

kable(filter_for_value) %>%
  kable_styling() %>%
  scroll_box(width = "700px", height = "400px")

#filter data set based on values that do not match the specified pattern (by using the minus symbol)
filter_for_not_a_value <- growth_data[-grepl("Strain", growth_data$Condition),]

kable(filter_for_not_a_value) %>%
  kable_styling() %>%
  scroll_box(width = "700px", height = "400px")

Other comments:

• Using regular expressions (programming symbol pattern) will increase their functionality
• Specified patterns are case sensitive (“t” does not equal “T”)
• Any matching pattern will be returned despite the context in which that pattern is located (i.e., grep(“the”, data) with return matches for “the“, “theme”, “heather”, “breathe“, and so on–this is where regular expressions are useful for specifying where in a string the pattern should appear.

Regular expressions are explained in the AWK and General notebook. But just to give an example lets just grep Strains that have a 3 letter number

filter_3letters <- growth_data[grepl("Strain[0-9][0-9][0-9]", growth_data$Condition),]

#check the structure of our data
kable(filter_3letters) %>%
  kable_styling() %>%
  scroll_box(width = "700px", height = "400px")

Here, [0-9] searches for every number from 0-9 and we look for three numbers. One important thing with subsetting is that the levels are still kept. So with the command above we remove every row that is not Strain101 and Strain230. However, the levels are still kept. Let’s check this:

filter_3letters$Condition

 [1] Strain101 Strain101 Strain101 Strain101 Strain101 Strain101 Strain101
 [8] Strain101 Strain101 Strain101 Strain101 Strain101 Strain101 Strain230
[15] Strain230 Strain230 Strain230 Strain230 Strain230 Strain230 Strain230
[22] Strain230 Strain230 Strain230 Strain101 Strain101 Strain101 Strain101
[29] Strain101 Strain101 Strain101 Strain101 Strain101 Strain101 Strain101
[36] Strain101 Strain230 Strain230 Strain230 Strain230 Strain230 Strain230
[43] Strain230 Strain230 Strain230 Strain230 Strain230 Strain230 Strain230
Levels: MgCl Strain101 Strain230 Strain28

So we see that MgCl and Strain28 are still in the levels even if they do not occur in our table itself. Sometimes when working with subsetted dataframes, i.e. when doing stats or plotting, this can interfer with our analysis. Here, it is useful to drop empty levels.

filter_3letters_clean <- droplevels(filter_3letters)
filter_3letters_clean$Condition

 [1] Strain101 Strain101 Strain101 Strain101 Strain101 Strain101 Strain101
 [8] Strain101 Strain101 Strain101 Strain101 Strain101 Strain101 Strain230
[15] Strain230 Strain230 Strain230 Strain230 Strain230 Strain230 Strain230
[22] Strain230 Strain230 Strain230 Strain101 Strain101 Strain101 Strain101
[29] Strain101 Strain101 Strain101 Strain101 Strain101 Strain101 Strain101
[36] Strain101 Strain230 Strain230 Strain230 Strain230 Strain230 Strain230
[43] Strain230 Strain230 Strain230 Strain230 Strain230 Strain230 Strain230
Levels: Strain101 Strain230

To do some basic stats we use the nutrient growth data set, which you find in the R_exercises folder.

• means() = print the mean of a certain column
• summary()= print the mean, median, max, min, … of a certain column.

#read in our data
growth_data <- read.table("Input_files/Growth_Data.txt", sep="\t", header=T,  quote = "")

#get the menas of our root length
mean(growth_data$Rootlength_cm)

[1] 9.046117

#return summary stats for the root length (using two different ways to write the syntax)
summary(growth_data$Rootlength_cm)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  4.236   6.721   8.099   9.046  11.793  14.716 

We can return the summary stats also for more than one column but for this we have to work with our square brackets again.

#return summary stats for the root length and shoot fresh weight
summary(growth_data[, c("Rootlength_cm", "FW_shoot_mg")])

 Rootlength_cm     FW_shoot_mg   
 Min.   : 4.236   Min.   : 3.75  
 1st Qu.: 6.721   1st Qu.: 7.69  
 Median : 8.099   Median :11.62  
 Mean   : 9.046   Mean   :15.23  
 3rd Qu.:11.793   3rd Qu.:22.05  
 Max.   :14.716   Max.   :39.39  

The function table builds a contingency table of the counts at each factor level, or combinations of factor levels. This is quite useful to count the number of data points with certain metadata associated. So for example for our genome data we can ask how many measurements we made for each condition.

#summarize how many measurements we have for each treatment
table(growth_data$SampleID)

    noP noP_101 noP_230  noP_28       P   P_101   P_230    P_28 
     13      13      11      15      15      12      13      13 

Here, we see that we have a slightly different number of measurements for each condition. The differences are slight, so we can ignore them, but this might be relevant if we have huge differences for some statistical analyses. This approach allows you to for larger datasets to easily check for samples that might be outliers in terms of measurements per sample.

The table command can get rather slow and there are some useful packages to speed up things and run more complicated mathematical operations. One example is the ddply() function of the plyr package. A useful feature of ddply is that this tool stores the output as a dataframe instead of a table.

Below we see examples were we summarize the data for the root length across different nutrient conditions. I.e. we do not want to have the mean for all root lengths but we want to see if roots have different lengths under low and high P conditions.

#load our package
library(plyr)

#calculate the mean root length across our Sample IDs
growth_data_summarized <- ddply(growth_data, .(Nutrient), summarize, RootLength = mean(Rootlength_cm))

#view data
kable(head(growth_data_summarized), format='markdown')

The structure of the command is as follows:

ddply(Input_Data, .(Colum_we_want_to_sum_arcoss), summarize, New_Column_name = mean(Column_we_want_to_calc_the_mean_for))

We can also summarize the data combining different conditions (i.e. nutrient and condition).

#Summarize our data across Nutrients and Conditions
growth_data_summarized <- ddply(growth_data, .(Nutrient, Condition), summarize, RootLength = mean(Rootlength_cm))

#view data
kable(growth_data_summarized, format='markdown')

We can also summarize the data for both root length and shoot weight.

#we can summarize more than one thing in one line of code
growth_data_summarized <- ddply(growth_data, .(Nutrient, Condition), summarize, RootLength = mean(Rootlength_cm), Freshweight = mean(FW_shoot_mg))

#view data
kable(growth_data_summarized, format='markdown')

We can also calculate the mean, sd and se in one line of code

#and we can use even fancier functions (i.e. get the se and sd), check the plyr package for more details
growth_data_summarized <- ddply(growth_data, .(Nutrient, Condition), summarize, RootLength = mean(Rootlength_cm), sd = sd (Rootlength_cm), se = sd(Rootlength_cm) / sqrt((length(Rootlength_cm))))

#view data
kable(growth_data_summarized, format='markdown')

unique allows to determine duplicate rows and allows us to subset our data for certain categories. For example, for very large dataframes we often can simplify things.

Here, if we have a lot of treatments and did the experiment a long time ago, we might want to ask for a table that lists the treatments.

#make unique contig list that still contains info of our bin ID
mapping_file <- unique(growth_data[,c("SampleID", "Nutrient", "Condition")])

#view data
kable((mapping_file), format='markdown')

We can also add additional metadata to our growth data.

One way to do this is the cbind() or rbind() functions. However, these functions require the two dataframes to have the exact number of columns or rows, which we do not have.

Here, the merge() function of the data.table package is very useful to merge data with different dimensions as long as they have a common pattern (i.e. the SampleID).

First lets build an artificial mapping file that incldues the number of days we grew our plants:

#make mapping that contains our basic sample info
mapping_file <- unique(growth_data[,c("SampleID", "Nutrient", "Condition")])

#add a new column, where we list the number of days the experiment ran
mapping_file$Days <- 12

#view data
kable(head(mapping_file), format='markdown')

Now we can use this mapping file and merge it with our growth data as follows:

#load our package
library(plyr)

#merge our mapping file with our growth data
new_data_frame <- merge(growth_data, mapping_file, by = "SampleID")

#view data
kable(head(new_data_frame), format='markdown')

This now is a good example to check that all went fine and that the new dataframes has the same number of rows (=measurements) compared to the original dataframe.

#control that all went fine
dim(growth_data)

[1] 105   5

dim(new_data_frame)

[1] 105   8

#if there is no match between dataframe 1 and dataframe 2 columns will by default be deleted. If you want to keep all columns do:
#new_data_frame <- merge(growth_data, mapping_file, by = "SampleID". all.x = T)

With dim we see that we still have 105 rows (i.e. measurements) and that we now added 3 new columns.

While this gets more difficult to read, sometimes it might be useful to combine several commands into one go to condense code Generally, it is easier to just write line by line especially if you read your code months later.

What we want to do: - in the example above we duplicate the columns for Nutrient and Condition and before merging we might first subset the mapping file to only include the info we want to merge. - so our two steps are: a) trim the mapping file and b) merge

To do this, we use these two lines of code:

#make mapping file more simple
mapping_reduced <- mapping_file[,c("SampleID", "Days")]

#merge
new_data_frame <- merge(growth_data, mapping_reduced, by = "SampleID")
kable(head(new_data_frame), format='markdown')

Now, this worked fine but requires a bit more code and we need to create one more object.

We could also combine these two kines of code into one line by subsetting our mapping file INSIDE the merge function as follows:

new_data_frame <- merge(growth_data, mapping_file[,c("SampleID", "Days")], by = "SampleID")
kable(head(new_data_frame), format='markdown')

We can also add new columns into original table, i.e. if we want to not show the fresh weight in mg but in g.

Below you can see that if we have numerical data in a column, we can use normal math operators (like +,-,/)

#convert mg to g and make a new column
growth_data$FW_shoot_g <- growth_data$FW_shoot_mg/10
kable(head(growth_data), format='markdown')

#we can also round our data
growth_data$FW_shoot_g <- round(growth_data$FW_shoot_mg/10, digits = 2)
kable(head(growth_data), format='markdown')

We can also add (or substract, etc …) values from different columns. I.e. here we could calculate the ratio from shoot weight to fresh weight.

#we can also do math with the values in two columns, i.e. if we want to calculate the ration between root length and fresh weight
growth_data$ratio <- growth_data$FW_shoot_mg/growth_data$Rootlength_cm
kable(head(growth_data), format='markdown')

Tapply is a base R function and allows to apply a function to each cell of a ragged array, that is to each (non-empty) group of values given by a unique combination of the levels of certain factors. It can be used as an alternative to ddply.

Suppose now we want to estimate the mean root length for each growth condition. Notice, how we can vary the index?

#index = factors we want to select
tapply(X = growth_data$Rootlength_cm, INDEX = growth_data$Nutrient, FUN = mean,na.rm = TRUE)

      noP         P 
 6.530845 11.513931 

#same but for two-way tables (this is not so useful here, but might be handy when you have different conditions for the same organism or a timecourse)
tapply(growth_data$Rootlength_cm, INDEX = list(growth_data$Nutrient, growth_data$Condition),FUN = mean, na.rm = TRUE)

        MgCl Strain101 Strain230  Strain28
noP 5.805198  7.094616  7.172883  6.200309
P   9.853457 12.729133 11.814678 12.007390

It is very simple to investigate linear relationships amongst variables in R. We want to estimate how a quantitative dependent variable changes according to the levels of one or more categorical independent variables. In the command below, the linear relationship between Rootlength_cm (the dependent variable, i.e. the one we’re trying to predict) and FW_shoot_mg (the independent variable or the predictor) is calculated.

We need to use the function summary to see the results of that command; coef extracts the best-fit coefficients, anova performs an analysis of variance; there are many other extractor functions.

#is there a correlation between freshweight and root length?
linearmodel <- lm(Rootlength_cm ~ FW_shoot_mg, data = growth_data)
linearmodel

Call:
lm(formula = Rootlength_cm ~ FW_shoot_mg, data = growth_data)

Coefficients:
(Intercept)  FW_shoot_mg  
     5.2912       0.2465  

#let's extract the entire t-table
summary(linearmodel) 

Call:
lm(formula = Rootlength_cm ~ FW_shoot_mg, data = growth_data)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.5263 -0.9978 -0.0181  1.0094  3.9292 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.29124    0.30469   17.37   <2e-16 ***
FW_shoot_mg  0.24648    0.01712   14.40   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.615 on 103 degrees of freedom
Multiple R-squared:  0.6681,    Adjusted R-squared:  0.6649 
F-statistic: 207.4 on 1 and 103 DF,  p-value: < 2.2e-16

Here one of the values is the model p-Value (bottom last line) and the p-Value of individual predictor variables (extreme right column under ‘Coefficients’). The p-Values are very important because, we can consider a linear model to be statistically significant only when both these p-Values are less that the pre-determined statistical significance level, which is ideally 0.05. This is visually interpreted by the significance stars at the end of the row. The more the stars beside the variable’s p-Value, the more significant the variable.

• Residuals: The section summarizes the residuals, the error between the prediction of the model and the actual results. Smaller residuals are better.
• Coefficients: For each variable and the intercept, a weight is produced and that weight has other attributes like the standard error, a t-test value and significance.
• Estimate: This is the weight given to the variable. In the simple regression case (one variable plus the intercept), for every increase in root length, the model predicts an increase of 0.24.
• Std. Error: Tells you how precisely was the estimate measured. It’s really only useful for calculating the t-value.
• t-value and Pr(>[t]): The t-value is calculated by taking the coefficient divided by the Std. Error. It is then used to test whether or not the coefficient is significantly different from zero. If it isn’t significant, then the coefficient really isn’t adding anything to the model and could be dropped or investigated further. Pr(>|t|) is the significance level.

Performance Measures:

Three sets of measurements are provided.

• Residual Standard Error: This is the standard deviation of the residuals. Smaller is better.
• Multiple / Adjusted R-Square: For one variable, the distinction doesn’t really matter. R-squared shows the amount of variance explained by the model. Adjusted R-Square takes into account the number of variables and is most useful for multiple-regression.
• F-Statistic: The F-test checks if at least one variable’s weight is significantly different than zero. This is a global test to help asses a model. If the p-value is not significant (e.g. greater than 0.05) than your model is essentially not doing anything.

We also can just print parts of the data:

#print only the coefficients 
coef(summary(linearmodel))

             Estimate Std. Error  t value     Pr(>|t|)
(Intercept) 5.2912436 0.30469198 17.36588 2.278273e-32
FW_shoot_mg 0.2464783 0.01711628 14.40023 2.039703e-26

#print only the anova stats
anova(linearmodel)

Analysis of Variance Table

Response: Rootlength_cm
             Df Sum Sq Mean Sq F value    Pr(>F)    
FW_shoot_mg   1 540.99  540.99  207.37 < 2.2e-16 ***
Residuals   103 268.71    2.61                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#plot add the best line to a plot
with(growth_data, plot(Rootlength_cm ~ FW_shoot_mg, col = 2))
abline(linearmodel)

If we look at the stats and the p value we see a nice correlation but also that we have two distinct clusters as well as more spread in the cluster that is more to the right. These clusters likely are the two different nutrient conditions and sometimes it might make sense to separate data to get a clearer picture. Something else to consider is to ask whether the data is normaly distributed and based on that what statistical test to choose.

Anyone can fit a linear model in R. The real test is analyzing the residuals (the error or the difference between actual and predicted results).

There are four things we’re looking for when analyzing residuals.

• The mean of the errors is zero (and the sum of the errors is zero)
• The distribution of the errors are normal.
• All of the errors are independent.
• Variance of errors is constant (Homoscedastic)

In R, you pull out the residuals by referencing the model and then the resid variable inside the model. Using the simple linear regression model (simple.fit) we’ll plot a few graphs to help illustrate any problems with the model.

Below some examples:

simple.fit <- linearmodel

layout(matrix(c(1,1,2,3),2,2,byrow=T))

#Rootlength_cm x Residuals Plot
plot(simple.fit$resid~growth_data$Rootlength_cm[order(growth_data$Rootlength_cm)],
 main="Rootlength_cm x Residuals\nfor Simple Regression",
 xlab="Marketing Rootlength_cm", ylab="Residuals")
abline(h=0,lty=2)

#Histogram of Residuals
hist(simple.fit$resid, main="Histogram of Residuals",
 ylab="Residuals")

#Q-Q Plot
qqnorm(simple.fit$resid)
qqline(simple.fit$resid)

The histogram and QQ-plot are the ways to visually evaluate if the residual fit a normal distribution.

• If the histogram looks like a bell-curve it might be normally distributed.
• If the QQ-plot has the vast majority of points on or very near the line, the residuals may be normally distributed.

Visualize our data via density plots

There are different ways to visualize this, one example is ggdensity of the ggpubr package.

library("ggpubr")

#is the data for my different variables normaly distributed
ggdensity(growth_data$Rootlength_cm)

We see nicely that we have two tails that likely represent the two nutrient conditions. To test this, we can simply subset the data as we have done before.

#subset
noP_data <- growth_data[growth_data$Nutrient=="noP",]

#plot
ggdensity(noP_data$Rootlength_cm)

No we see that indeed the tool tails we see are seem to be due to our two nutrient conditions.

Visualize our data via Q-Q plots

Another way to represent data is in a Q-Q plot: Q-Q plot (or quantile-quantile plot) draws the correlation between a given sample and the normal distribution. A 45-degree reference line is also plotted.

#plot all data
ggqqplot(growth_data$Rootlength_cm)

#only plot low P data
ggqqplot(noP_data$Rootlength_cm)

Again, here we see that our data for the indivdual growth conditions fit quite nicely into normal distribution.

Test for normality

#for all data
shapiro.test(growth_data$Rootlength_cm)

    Shapiro-Wilk normality test

data:  growth_data$Rootlength_cm
W = 0.91773, p-value = 6.744e-06

#for noP only
shapiro.test(noP_data$Rootlength_cm)

    Shapiro-Wilk normality test

data:  noP_data$Rootlength_cm
W = 0.96403, p-value = 0.1171

The shapiro.test tests the NULL hypothesis that the samples came from a Normal distribution. This means that if your p-value <= 0.05, then you would reject the NULL hypothesis that the samples came from a normal distribution.

From the output, the p-value < 0.05 for our complete dataset implies that the distribution of the data is significantly different from normal distribution. In other words, we can not assume the normality. However, we expect quite some differences dependent on the growth pattern and once we only look at our low P data we see that our indeed is normally distributed.

Notice: Shapiro works only for sample sizes between 3-5000 numbers since when you feed it more data, the chances of the null hypothesis being rejected becomes larger. An alternative is the Anderson-Darling test that however has a similar problem with the Shapiro Wilk test. For large samples, you are most likely to reject the null hypothesis, so be aware of this.

library(nortest)
ad.test(noP_data$Rootlength_cm)

    Anderson-Darling normality test

data:  noP_data$Rootlength_cm
A = 0.58838, p-value = 0.1187

When we visuallize the data, we see that there is a difference between the nutrient conditions but we want to know whether it is significant and more importantly, whether there is also a difference based on our treatments with different strains of microbes.

#lets visually compare our data with ggpubr again
library("ggpubr")

ggboxplot(growth_data, x = "Condition", y = "Rootlength_cm", color = "Nutrient",
          palette = c("#00AFBB", "#E7B800"))

# We want to know whether root length depends on nutrient treatment
aov <- aov(Rootlength_cm ~ Nutrient, data = growth_data)
summary(aov)

             Df Sum Sq Mean Sq F value Pr(>F)    
Nutrient      1  651.8   651.8     425 <2e-16 ***
Residuals   103  157.9     1.5                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Here, we see that there are significant differences based on our nutrient treatments. Now lets see how we can look at both the nutrient treatment and growth conditions.

# We want to know if root length depends on condition and nutrient
aov <- aov(Rootlength_cm ~ Nutrient + Condition, data = growth_data)
summary(aov)

             Df Sum Sq Mean Sq F value   Pr(>F)    
Nutrient      1  651.8   651.8  712.58  < 2e-16 ***
Condition     3   66.5    22.2   24.23 7.27e-12 ***
Residuals   100   91.5     0.9                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

From the ANOVA table we can conclude that both nutrient conditon and treatment are statistically significant. Nutrient treatment is the most significant factor variable.

Not the above fitted model is called additive model. It makes an assumption that the two factor variables are independent. If you think that these two variables might interact to create an synergistic effect, replace the plus symbol (+) by an asterisk (*), as follow.

# Two-way ANOVA with interaction effect

# These two calls are equivalent
aov <- aov(Rootlength_cm ~ Nutrient * Condition, data = growth_data)
aov <- aov(Rootlength_cm ~ Nutrient + Condition + Nutrient:Condition, data = growth_data)

#summarize the aov
summary(aov)

                   Df Sum Sq Mean Sq F value   Pr(>F)    
Nutrient            1  651.8   651.8 817.131  < 2e-16 ***
Condition           3   66.5    22.2  27.780 4.72e-13 ***
Nutrient:Condition  3   14.1     4.7   5.891 0.000976 ***
Residuals          97   77.4     0.8                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

It can be seen that the two main effects (supp and dose) are statistically significant, as well as their interaction.

Note that, in the situation where the interaction is not significant you should use the additive model.

In ANOVA test, a significant p-value indicates that some of the group means are different, but we don’t know which pairs of groups are different. It’s possible to perform multiple pairwise-comparison, to determine if the mean difference between specific pairs of group are statistically significant.

As the ANOVA test is significant, we can compute Tukey HSD (Tukey Honest Significant Differences). Tukey test is a single-step multiple comparison procedure and statistical test. It is a post-hoc analysis, what means that it is used in conjunction with an ANOVA.

#test with anova
aov <- aov(Rootlength_cm ~ Nutrient * Condition, data = growth_data)

#run tukey
TukeyHSD(aov)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Rootlength_cm ~ Nutrient * Condition, data = growth_data)

$Nutrient
          diff      lwr      upr p adj
P-noP 4.983086 4.637105 5.329067     0

$Condition
                          diff        lwr        upr     p adj
Strain101-MgCl       2.1029052  1.4604930  2.7453173 0.0000000
Strain230-MgCl       1.6836191  1.0341749  2.3330633 0.0000000
Strain28-MgCl        1.2784796  0.6545138  1.9024455 0.0000034
Strain230-Strain101 -0.4192861 -1.0864726  0.2479004 0.3597470
Strain28-Strain101  -0.8244256 -1.4668377 -0.1820134 0.0061445
Strain28-Strain230  -0.4051395 -1.0545837  0.2443048 0.3663501

$`Nutrient:Condition`
                                  diff        lwr        upr     p adj
P:MgCl-noP:MgCl              4.0482585  2.9997859  5.0967311 0.0000000
noP:Strain101-noP:MgCl       1.2894180  0.2041458  2.3746903 0.0087961
P:Strain101-noP:MgCl         6.9239348  5.8162834  8.0315861 0.0000000
noP:Strain230-noP:MgCl       1.3676850  0.2341551  2.5012149 0.0073003
P:Strain230-noP:MgCl         6.0094793  4.9242070  7.0947515 0.0000000
noP:Strain28-noP:MgCl        0.3951101 -0.6533625  1.4435827 0.9391499
P:Strain28-noP:MgCl          6.2021911  5.1169189  7.2874633 0.0000000
noP:Strain101-P:MgCl        -2.7588405 -3.8073131 -1.7103679 0.0000000
P:Strain101-P:MgCl           2.8756763  1.8040558  3.9472967 0.0000000
noP:Strain230-P:MgCl        -2.6805735 -3.7789218 -1.5822252 0.0000000
P:Strain230-P:MgCl           1.9612207  0.9127482  3.0096933 0.0000023
noP:Strain28-P:MgCl         -3.6531484 -4.6634819 -2.6428149 0.0000000
P:Strain28-P:MgCl            2.1539326  1.1054600  3.2024052 0.0000002
P:Strain101-noP:Strain101    5.6345168  4.5268654  6.7421681 0.0000000
noP:Strain230-noP:Strain101  0.0782670 -1.0552629  1.2117969 0.9999989
P:Strain230-noP:Strain101    4.7200613  3.6347890  5.8053335 0.0000000
noP:Strain28-noP:Strain101  -0.8943079 -1.9427805  0.1541647 0.1537008
P:Strain28-noP:Strain101     4.9127731  3.8275008  5.9980453 0.0000000
noP:Strain230-P:Strain101   -5.5562498 -6.7112241 -4.4012755 0.0000000
P:Strain230-P:Strain101     -0.9144555 -2.0221069  0.1931958 0.1845494
noP:Strain28-P:Strain101    -6.5288247 -7.6004451 -5.4572042 0.0000000
P:Strain28-P:Strain101      -0.7217437 -1.8293950  0.3859077 0.4749256
P:Strain230-noP:Strain230    4.6417943  3.5082644  5.7753242 0.0000000
noP:Strain28-noP:Strain230  -0.9725749 -2.0709232  0.1257734 0.1223122
P:Strain28-noP:Strain230     4.8345061  3.7009762  5.9680360 0.0000000
noP:Strain28-P:Strain230    -5.6143692 -6.6628417 -4.5658966 0.0000000
P:Strain28-P:Strain230       0.1927118 -0.8925604  1.2779841 0.9993212
P:Strain28-noP:Strain28      5.8070810  4.7586084  6.8555536 0.0000000

We can see that most differences are significant, with the exception of Strain28, which in most cases does not show an effect.

For some representations it is useful to plot significant letters. We can do this using some extra packages as follows:

#load library
library(agricolae)

#separate nutrient conditions
noP_data <- growth_data[growth_data$Nutrient == "noP", ]

#run an anova
aov_noP <- aov(Rootlength_cm ~ Condition, data = noP_data)

#run test
HSD.test(aov_noP,"Condition", group=TRUE,console=TRUE)

Study: aov_noP ~ "Condition"

HSD Test for Rootlength_cm 

Mean Square Error:  0.606883 

Condition,  means

          Rootlength_cm       std  r      Min      Max
MgCl           5.805198 0.7738124 13 4.236348 6.834720
Strain101      7.094616 0.5939804 13 5.947965 7.973207
Strain230      7.172883 0.6117203 11 6.120089 8.197116
Strain28       6.200309 0.9988989 15 4.384756 7.428674

Alpha: 0.05 ; DF Error: 48 
Critical Value of Studentized Range: 3.763749 

Groups according to probability of means differences and alpha level( 0.05 )

Treatments with the same letter are not significantly different.

          Rootlength_cm groups
Strain230      7.172883      a
Strain101      7.094616      a
Strain28       6.200309      b
MgCl           5.805198      b

Sometimes it is necessary to repeat a calculation multiple times, e.g. calculate the sum of each row of a matrix. You can use for loops to do this.

In the first loop, we:

• define that we want to loop through rows 1 to 5 (1:5)
• for (i in 1:5) –> we say that we want to run a for loop, in which we name our variable i
• {} –> defines the operation we want to loop through
• in our case we want to sum a specific row (rows 1 trhough 5)
• so in i we store the numbers 1,2,3,4 and 5 and then run through the sum functions 5x

#build a dataframe
m <- matrix(1:15, 5)
m

     [,1] [,2] [,3]
[1,]    1    6   11
[2,]    2    7   12
[3,]    3    8   13
[4,]    4    9   14
[5,]    5   10   15

#sum rows without a loop --> sum row 1
sum(m[1, ])

[1] 18

#sum rows without a loop --> sum row 2, ...
sum(m[2, ])

[1] 21

#run a loop to get the sum for all five rows
for (i in 1:5) {
print(sum(m[i, ]))
}

[1] 18
[1] 21
[1] 24
[1] 27
[1] 30

In the second loop we store our results in a new dataframe and do not just print them to the screen. Therefore, we:

• define a variable (named results) in which we store an empty vector
• we need this empty vector to have something to store our results in
• start the foor loop and define i as before
• run the sum as before but now store the results for each interation in results

#store the results, Note that here, the variable results had to be created before as an empty vector.
results <- c()
for (i in 1:5) {
results[i] <- sum(m[i, ])
}

results

[1] 18 21 24 27 30

Very often, there are some R packages or functions that are faster than loops. The apply() function is such an example. apply() applies a funcion, i.e. sum, across a matrix.

However, since loops are also a very useful feature in bash or python it is useful to understand their general concept.

#way to do the same using the apply function (faster than loops)
results2 <- apply(m, 1, sum)
results2

[1] 18 21 24 27 30

The first argument of apply() has to be a matrix, the second argument specifies the dimension of the matrix (1 means rows and 2 means columns) on which the function should be applied, and the third argument has to be a function that takes a vector as argument. Here, the function sum is applied to each row of m.

Actually, there is a family of apply functions, depending on the object you pass as an input, and/or the object you want as output. You can read a brief tutorial under this link: (http://datascienceplus.com/using-the-apply-family-of-functions-in-r/).

I.e. apply(data, 1, mean) = apply mean on each row apply(data, 2, mean) = apply mean on each column

Two other examples are sapply and lapply, which work on lists.

One can also create a decision making structure in R. A decision making structure has at least one condition to be evaluated together with a statement or statements to be evaluated if the condition is TRUE, and optionally, other statements to be executed if the condition is FALSE, as we can see in the figure.

The test_expression has to be logical, i.e., it has to be a expression that, when evaluated, is either TRUE or FALSE. The logical operators listed above can be used to construct them. For example, we can use an if-else statement to check if a number is positive or negative,

#store a random number
x <- -5

#write a loop and ask if our number is positive or negative
if (x > 0) {
print("Positive number")
} else if (x == 0) {
print("Number is zero")
} else {
print("Negative number")
}

[1] "Negative number"

Note that in the example there was an else if. In that way we can check more than 2 conditions.

This of course is a very simplistic example. But if-else statements can be useful if we want to deal with larger datatables, i.e. our growth data.

#apply ifelse
growth_data$category <- ifelse(growth_data$Rootlength_cm<=10.0, "small", "large")

#check the structure of our data
kable(growth_data) %>%
  kable_styling() %>%
  scroll_box(width = "700px", height = "400px")

Here, we:

• Create a new column, with the name category
• If a value in our Rootlength_cm column is equal to or smaller than 10, we want to say this category is small.
• If that is not the case, we want to say it is large

So the function works like this:

ifelse(our_test, value_to_return_if_test_is_true, value_to_return_if_test_is_false)

We can also combine diffent statements with the & (AND) or | (OR) symbol. I.e. we only send things in the big category if the roots are longer than 10cm AND the shoot weight is larger than 15mg

#apply ifelse
growth_data$category <- ifelse(growth_data$Rootlength_cm>10 & growth_data$FW_shoot_mg>15, "large", "small")

#check the structure of our data
kable(growth_data) %>%
  kable_styling() %>%
  scroll_box(width = "700px", height = "400px")